#This files takes functions from MCES..._functions and adapts to Mixed Logit case 
#when relevant.
#This version intended for "macro" estimate of eta

####################################
#Equilibrium p and s on each market#
####################################
eqbm.mlog <- function(dat,market,v=1,tar.change=0,verbose=FALSE) { 
  # the new version of nov 2020 can deal with multiple X and outside good
  # most important, it sources the main subs-functions. 
  #
  source("eqbm_MLOG_subfunctions.R",local=TRUE)
  #
  ###########################
  #start of eqbm computation#  
  df.m <- data.frame(subset(dat$models,n==market),tar.change=tar.change)
  df.m <- within(df.m,{tar.change[i==n]<-0})
  i <- df.m$i
  if(U0==0) xvars <- paste0("x",1:n.x) else xvars <- paste0("x",0:n.x)
  x <- as.matrix(df.m[,xvars]) # if df.m is only a data.frame
  colnames(x) <- xvars  # no effect for n.x>1, but names it x1 instead of just x
  xi <- df.m$xi
  mc <- df.m$mc * (df.m$tau + df.m$tar.change*df.m$tariff.change) # apply trade costs and name it mc, for consistency with subordinate functions
  #
  df <- subset(dat$households,n==market)
  if(U0==0) bvec <- paste0("b",1:n.x,".h") else bvec <- paste0("b",0:n.x,".h")
  b.h <- as.matrix(df[,bvec]) 
  colnames(b.h) <- bvec  # no effect for n.x>1, but names it b1.h instead of just b.h
  alpha.h <- df$alpha.h
  y.h <- df$y.h
  Y <- sum(df$y.h)
  co.own <- co.own.mat
  # FPI to get the MCES predict prices and market shares
  FPIresults <- fpiter(mc,fixptfn=price.update,nu=settings.var$speed[v],control=list(trace=verbose,maxiter=MAXIT))
  n.loops  <- FPIresults$fpevals
  p <- FPIresults$par
  s <- colMeans(prob.mat(p))
  pelas.m <- ownelas(prm=prob.mat(p),share=s,price=p) # model level own price elasticities
  superelas.m <- superelas(prm=prob.mat(p),share=s,price=p) # model level own super elasticity
  ptr.sp.m <- (pelas.m)/(pelas.m + 1 - superelas.m)
  pte.sp.m <- (pelas.m + 1)/(pelas.m + 1 - superelas.m)
  #average income per model and ISE
  avg.inc.m <- avg.inc(prob.mat(p))   
  avg.inc.results <- lm(log(avg.inc.m)~log(p))
  ISE <- as.numeric(avg.inc.results$coefficients[2])
  return(data.frame(m=df.m$m,f=df.m$f,i,n=market,x,xi,mc=df.m$mc,cit=mc,s,p,pelas.m,ptr.sp.m,pte.sp.m,avg.inc.m,ISE,n.loops=n.loops,Y=Y))
  }

#Putting all markets together
DGP.endog.mlog <- function(parallel=TRUE,v=1,tar.change=0,data.exog){ #get equilibrium market shares and prices for each market and append them
  if(parallel==T)  foreach(i=1:n.countries, .combine = 'rbind')  %dorng%  eqbm.mlog(dat=data.exog,i,v,tar.change) else foreach(i=1:n.countries, .combine = 'rbind')  %do%  eqbm.mlog(dat=data.exog,i,v,tar.change) 
}

##########################
######## MC part #########

#Monte Carlo repetitions MM - MC
monte.func.mlog.mm <- function(repi,settings,v=1,market=1) { # monte carlo stage
  data.exog.mm <- DGP.exog(settings=settings,v=v,grid=subset(tau.grid.allv,vtau==v,select=-vtau)) #gen exog data
  #pre
  eqbm.pre.mm  <- as.data.table(DGP.endog.mlog(parallel=F,v=v,tar.change=0,data.exog = data.exog.mm))
  eqbm.pre.mm.with.ces <- as.data.table(estapprox.mm(eqbm.pre.mm)) #add ces approximation results
  eqbm.pre.mm.with.ces <- as.data.table(estapproxagg.mm(eqbm.pre.mm.with.ces)) #add ces approximation results
  #post:
  eqbm.post <- eqbm.mlog(dat=data.exog.mm,market=market,v=v,tar.change=1)
  eqbm.pre.with.ces <- eqbm.pre.mm.with.ces[n==market,] #keep only one destination
  eqbm.pre <- eqbm.pre.mm[n==market,] #keep only one destination
  #AFTER THIS, nothing in the code is changed compared to 2ctry, since we have only one destination left.
  eqbm.post.with.ces <- as.data.table(CF.approx(eqbm.pre.with.ces,eqbm.post)) #add ces approximation results
  return(data.frame(repi=repi,v=v,CF.compare(eqbm.pre.with.ces,eqbm.post.with.ces)))
}

#Doing multiple monte carlo reps for a single market
monte.rep.mlog.mm <- function(parallel=TRUE,settings,v=1,market=1){ #get equilibrium market shares and prices for each market and append them
  if(parallel==T)  foreach(i=1:n.reps, .combine = 'rbind')  %dorng%  monte.func.mlog.mm(i,settings,v,market) else foreach(i=1:n.reps, .combine = 'rbind')  %do%  monte.func.mlog.mm(i,settings,v,market) 
}

